// Looking at historical sorting patterns of the US wrt time zones
//
// Jeff Shrader & Matthew Gibson
// Creation date: 2014-06-01
// Time-stamp: "2014-06-05 20:46:06 jgs"

// Preliminaries
clear

local work "/DIRECTORY"
cd `work'

capture log close
log using "`work'/logs/historical_sorting.log", replace
set more off

// Bring in census data
use "data/census_population.dta", clear
drop if meta_merge != 3
drop meta_merge
encode fips, gen(id)
tsset id year, delta(10)
bysort fips (year): gen pop_growth = (totpop - totpop[_n-1])/totpop[_n-1]

// Drawing rough time zone boundaries based on newspaper articles mentioned
// on statoids.com
gen time_zone_1883a = "E" if cnty_long < 84
replace time_zone_1883a = "C" if cnty_long < 102 & time_zone_1883a == ""
replace time_zone_1883a = "M" if cnty_long < 114 & time_zone_1883a == ""
replace time_zone_1883a = "P" if time_zone_1883a == ""
// Hawaii and Alaska are not states yet, so we don't need to worry about them!
gen tzdist_E_1883a = cnty_long - 67.28 if time_zone_1883a == "E"
replace tzdist_E_1883a = cnty_long - 84 if time_zone_1883a == "C"
replace tzdist_E_1883a = cnty_long - 102 if time_zone_1883a == "M"
replace tzdist_E_1883a = cnty_long - 114 if time_zone_1883a == "P"
gen tzdist_W_1883a = 84 - cnty_long  if time_zone_1883a == "E"
replace tzdist_W_1883a = 102 - cnty_long if time_zone_1883a == "C"
replace tzdist_W_1883a = 114 - cnty_long if time_zone_1883a == "M"
replace tzdist_W_1883a = 125 - cnty_long if time_zone_1883a == "P"


// Creating medians
gen east_side = (tzdist_E_1883a <= 1.3)
gen west_side = (tzdist_W_1883a <= 1.72)
gen central = (east_side == 0 & west_side == 0)
bysort year: egen med_growth_E_a = median(pop_growth) if east_side == 1 
bysort year: egen med_growth_W_a = median(pop_growth) if west_side == 1
bysort year: egen med_growth_C_a = median(pop_growth) if central == 1 


sort id year
gen d_growth = pop_growth - L.pop_growth


bysort id: egen growth_1880_1850 = mean(pop_growth) if year >= 1850 & year <= 1880
bysort id: egen growth_1920_1890 = mean(pop_growth) if year >= 1890 & year <= 1920
spread growth_1880_1850, by(id)
spread growth_1920_1890, by(id)
gen d_growth_long = growth_1880_1850 - growth_1920_1890


spread med_growth_*_a, by(year)
preserve
bysort year: gen first = _n
keep if first == 1
twoway (lpoly med_growth_E_a year if year <= 1880, bwidth(10) lcolor(red)) ///
  (lpoly med_growth_E_a year if year > 1880,  bwidth(10) lcolor(red)) ///
  (lpoly med_growth_W_a year if year <= 1880,  bwidth(10) lcolor(blue)) ///
  (lpoly med_growth_W_a year if year > 1880,  bwidth(10) lcolor(blue)) ///
  (scatter med_growth_E_a med_growth_W_a year, mcolor(red blue)) ///
  (pci 0 1883 .6 1883, lcolor(black)) if first == 1 & year < 1960, ///
  legend(label(1 "Eastern side of TZ") label(3 "Western side of TZ") ///
  order(1 3)) ytitle("Pop. growth") 
graph export "graphs/historical_sorting/median_diff_1883a.pdf", as(pdf) replace

twoway (line med_growth_E_a year if year <= 1918,  lcolor(red)) ///
  (line med_growth_E_a year if year > 1918,  lcolor(red)) ///
  (line med_growth_W_a year if year <= 1918, lcolor(blue)) ///
  (line med_growth_W_a year if year > 1918, lcolor(blue)) ///
  (scatter med_growth_E_a med_growth_W_a year, mcolor(red blue)) ///
  (pci 0 1918 .6 1918, lcolor(black)) if first == 1 & year <= 1950 & year >= 1800, ///
  legend(label(1 "Eastern side of 1918 time zone") label(3 "Western side of 1918 time zone") ///
  order(1 3)) ytitle("Population growth") 
graph save "graphs/historical_sorting/median_diff_1918.gph", replace
restore


preserve
bysort year: gen first = _n
keep if first == 1
twoway (lpoly med_growth_E_a year if year <= 1880, bwidth(10) lcolor(red)) ///
  (lpoly med_growth_E_a year if year > 1880 & year < 1920,  bwidth(10) lcolor(red)) ///
  (lpoly med_growth_E_a year if year >= 1920,  bwidth(10) lcolor(red)) ///
  (lpoly med_growth_W_a year if year <= 1880,  bwidth(10) lcolor(blue)) ///
  (lpoly med_growth_W_a year if year > 1880 & year < 1920,  bwidth(10) lcolor(blue)) ///
  (lpoly med_growth_W_a year if year > 1910,  bwidth(10) lcolor(blue)) ///
  (scatter med_growth_E_a med_growth_W_a year, mcolor(red blue)) ///
  (pci 0 1883 .6 1883, lcolor(black)) ///
  (pci 0 1918 .6 1918, lcolor(black)) if first == 1, ///
  legend(label(1 "Eastern side of TZ") label(3 "Western side of TZ") ///
  order(1 3)) ytitle("Pop. growth") 
graph export "graphs/historical_sorting/median_diff_both_1883a.pdf", as(pdf) replace
restore

// Alternative 1883 boundaries based on 15 degree increments
// Note that this is quite different than the above the the farthest
// east and west time zones, since we have fixed boundaries, even
// if they are in the middle of the ocean.
gen time_zone_1883b = "A" if cnty_long < 75
replace time_zone_1883b = "E" if cnty_long < 90 & cnty_long >= 75
replace time_zone_1883b = "C" if cnty_long < 105 & cnty_long >= 90
replace time_zone_1883b = "M" if cnty_long < 120 & cnty_long >= 105
replace time_zone_1883b = "P" if cnty_long >= 120
gen tzdist_E_1883b = cnty_long - 60 if time_zone_1883b == "A"
replace tzdist_E_1883b = cnty_long - 75 if time_zone_1883b == "E"
replace tzdist_E_1883b = cnty_long - 90 if time_zone_1883b == "C"
replace tzdist_E_1883b = cnty_long - 105 if time_zone_1883b == "M"
replace tzdist_E_1883b = cnty_long - 120 if time_zone_1883b == "P"
gen tzdist_W_1883b = 75 - cnty_long  if time_zone_1883b == "A"
replace tzdist_W_1883b = 90 - cnty_long if time_zone_1883b == "E"
replace tzdist_W_1883b = 105 - cnty_long if time_zone_1883b == "C"
replace tzdist_W_1883b = 120 - cnty_long if time_zone_1883b == "M"
replace tzdist_W_1883b = 135 - cnty_long if time_zone_1883b == "P"

// Creating medians
gen east_side_b = (tzdist_E_1883b <= 1.3)
gen west_side_b = (tzdist_W_1883b <= 1.8)
gen central_b = (east_side_b == 0 & west_side_b == 0)
bysort year: egen med_growth_E_b = median(pop_growth) if east_side_b == 1 & pop_growth < 5 & pop_growth > -5 
bysort year: egen med_growth_W_b = median(pop_growth) if west_side_b == 1 & pop_growth < 5 & pop_growth > -5 
bysort year: egen med_growth_C_b = median(pop_growth) if central_b == 1  & pop_growth < 5 & pop_growth > -5 

sort year


sort id year


// A better median graph
spread med_growth_*_b, by(year)
preserve
bysort year: gen first = _n
keep if first == 1
twoway (lpoly med_growth_E_b year if year <= 1880, bwidth(10) lcolor(red)) ///
  (lpoly med_growth_E_b year if year > 1880,  bwidth(10) lcolor(red)) ///
  (lpoly med_growth_W_b year if year <= 1880,  bwidth(10) lcolor(blue)) ///
  (lpoly med_growth_W_b year if year > 1880,  bwidth(10) lcolor(blue)) ///
  (scatter med_growth_E_b med_growth_W_b year, mcolor(red blue)) ///
  (pci 0 1883 .5 1883, lcolor(black)) if first == 1 & year < 1960, ///
  legend(label(1 "Eastern side of TZ") label(3 "Western side of TZ") ///
  order(1 3)) ytitle("Pop. growth") 
graph export "graphs/historical_sorting/median_diff_1883b.pdf", as(pdf) replace

twoway (line med_growth_E_b year if year <= 1880,  lcolor(red)) ///
  (line med_growth_E_b year if year > 1880,  lcolor(red)) ///
  (line med_growth_W_b year if year <= 1880, lcolor(blue)) ///
  (line med_growth_W_b year if year > 1880, lcolor(blue)) ///
  (scatter med_growth_E_b med_growth_W_b year, mcolor(red blue)) ///
  (pci 0 1883 .6 1883, lcolor(black)) if first == 1 & year < 1960 & year >= 1800, ///
  legend(label(1 "Eastern side of 1883 time zone") label(3 "Western side of 1883 time zone") ///
  order(1 3)) ytitle("Population growth") 
graph save "graphs/historical_sorting/median_diff_1883b.gph", replace
restore

preserve
bysort year: gen first = _n
keep if first == 1
twoway (lpoly med_growth_E_b year if year <= 1880, bwidth(10) lcolor(red)) ///
  (lpoly med_growth_E_b year if year > 1880 & year < 1920,  bwidth(10) lcolor(red)) ///
  (lpoly med_growth_E_b year if year >= 1920,  bwidth(10) lcolor(red)) ///
  (lpoly med_growth_W_b year if year <= 1880,  bwidth(10) lcolor(blue)) ///
  (lpoly med_growth_W_b year if year > 1880 & year < 1920,  bwidth(10) lcolor(blue)) ///
  (lpoly med_growth_W_b year if year > 1910,  bwidth(10) lcolor(blue)) ///
  (scatter med_growth_E_b med_growth_W_b year, mcolor(red blue)) ///
  (pci 0 1883 .6 1883, lcolor(black)) ///
  (pci 0 1918 .6 1918, lcolor(black)) if first == 1, ///
  legend(label(1 "Eastern side of TZ") label(3 "Western side of TZ") ///
  order(1 3)) ytitle("Pop. growth") 
graph export "graphs/historical_sorting/median_diff_both_1883b.pdf", as(pdf) replace
restore

// Simple t-test of change in growth around policy date
gen e_w = 1 if east_side_b == 1
replace e_w = 0 if west_side_b == 1
ttest d_growth if year == 1890, by(e_w)


// Handling the westward expansion in some way
reg pop_growth cnty_longitude
predict pg_res, resid
bysort year: egen med_gr_E_b = median(pg_res) if east_side_b == 1 
bysort year: egen med_gr_W_b = median(pg_res) if west_side_b == 1
bysort year: egen med_gr_C_b = median(pg_res) if central_b == 1 

// A better median graph
spread med_gr_*_b, by(year)
preserve
bysort year: gen first = _n
keep if first == 1
twoway (lpoly med_gr_E_b year if year <= 1880, bwidth(10) lcolor(red)) ///
  (lpoly med_gr_E_b year if year > 1880,  bwidth(10) lcolor(red)) ///
  (lpoly med_gr_W_b year if year <= 1880,  bwidth(10) lcolor(blue)) ///
  (lpoly med_gr_W_b year if year > 1880,  bwidth(10) lcolor(blue)) ///
  (scatter med_gr_E_b med_gr_W_b year, mcolor(red blue)) ///
  (pci 0 1883 .5 1883, lcolor(black)) if first == 1 & year < 1960, ///
  legend(label(1 "Eastern side of TZ") label(3 "Western side of TZ") ///
  order(1 3)) ytitle("Pop. growth") 
graph export "graphs/historical_sorting/median_diff_resid_1883b.pdf", as(pdf) replace
restore

preserve
bysort year: gen first = _n
keep if first == 1
twoway (lpoly med_gr_E_b year if year <= 1880, bwidth(10) lcolor(red)) ///
  (lpoly med_gr_E_b year if year > 1880 & year < 1920,  bwidth(10) lcolor(red)) ///
  (lpoly med_gr_E_b year if year >= 1920,  bwidth(10) lcolor(red)) ///
  (lpoly med_gr_W_b year if year <= 1880,  bwidth(10) lcolor(blue)) ///
  (lpoly med_gr_W_b year if year > 1880 & year < 1920,  bwidth(10) lcolor(blue)) ///
  (lpoly med_gr_W_b year if year > 1910,  bwidth(10) lcolor(blue)) ///
  (scatter med_gr_E_b med_gr_W_b year, mcolor(red blue)) ///
  (pci 0 1883 .6 1883, lcolor(black)) ///
  (pci 0 1918 .6 1918, lcolor(black)) if first == 1, ///
  legend(label(1 "Eastern side of TZ") label(3 "Western side of TZ") ///
  order(1 3)) ytitle("Pop. growth") 
graph export "graphs/historical_sorting/median_diff_resid_both_1883b.pdf", as(pdf) replace
restore



// Current boundaries
save Data/temp.dta, replace
use Data/temp.dta, clear
sort fips
merge m:1 fips using "data/county_tz/county_time_zone.dta"
drop if _merge == 2
drop _merge

rename time_zone time_zone_curr
gen tzdist_E_curr = cnty_long - 67.28 if time_zone_curr == "E"
replace tzdist_E_curr = cnty_long - 84.5 if time_zone_curr == "C"
replace tzdist_E_curr = cnty_long - 100.2 if time_zone_curr == "M"
replace tzdist_E_curr = cnty_long - 114.28 if upper(time_zone_curr) == "P"
gen tzdist_W_curr = 89.2 - cnty_long if time_zone_curr == "E"
replace tzdist_W_curr = 104.5 - cnty_long if time_zone_curr == "C"
replace tzdist_W_curr = 117.15 - cnty_long if time_zone_curr == "M"
replace tzdist_W_curr = 124.25 - cnty_long if upper(time_zone_curr) == "P"

// Creating medians
gen east_side_c = (tzdist_E_curr <= 3.6)
gen west_side_c = (tzdist_W_curr <= 4.1)
gen central_c = (east_side_c == 0 & west_side_c == 0)
bysort year: egen med_growth_E_c = median(pop_growth) if east_side_c == 1 & pop_growth < 5 & pop_growth > -5 
bysort year: egen med_growth_W_c = median(pop_growth) if west_side_c == 1 & pop_growth < 5 & pop_growth > -5 
bysort year: egen med_growth_C_c = median(pop_growth) if central_c == 1  & pop_growth < 5 & pop_growth > -5 

sort year
twoway (line med_growth_*_c year) (pci 0 1883 .6 1883), legend(label(1 "Eastern side of TZ") label(2 "Western side of TZ") label(3 "Center of TZ") order(1 2 3)) ytitle("Pop. growth") title("Median growth rates, 15deg")
graph export "graphs/historical_sorting/median_growth_curr.pdf", as(pdf) replace

sort id year
twoway (lpolyci d_growth tzdist_E_curr) if d_growth > -3 & d_growth < 2 & year == 1890, legend(off) xtitle("Distance from E. TZ boundary") ytitle("Change in pop. growth") title("Change in growth rate 1880-1890, 15deg")
graph export "graphs/historical_sorting/dgrowth_1890b.pdf", as(pdf) replace

twoway (lpolyci d_growth_long tzdist_E_curr) if d_growth > -3 & d_growth < 2 & year == 1890, legend(off) xtitle("Distance from E. TZ boundary") ytitle("Change in pop. growth") title("Change in average growth rate before and after 1883, 15deg")
graph export "graphs/historical_sorting/dgrowth_long_1890b.pdf", as(pdf) replace

// A better median graph
spread med_growth_*_c, by(year)
preserve
bysort year: gen first = _n
keep if first == 1
twoway (lpoly med_growth_E_c year if year <= 1880, bwidth(10) lcolor(red)) ///
  (lpoly med_growth_E_c year if year > 1880,  bwidth(10) lcolor(red)) ///
  (lpoly med_growth_W_c year if year <= 1880,  bwidth(10) lcolor(blue)) ///
  (lpoly med_growth_W_c year if year > 1880,  bwidth(10) lcolor(blue)) ///
  (scatter med_growth_E_c med_growth_W_c year, mcolor(red blue)) ///
  (pci 0 1883 .5 1883, lcolor(black)) if first == 1 & year < 1960, ///
  legend(label(1 "Eastern side of TZ") label(3 "Western side of TZ") ///
  order(1 3)) ytitle("Pop. growth") 
graph export "graphs/historical_sorting/median_diff_curr.pdf", as(pdf) replace
restore

preserve
bysort year: gen first = _n
keep if first == 1
twoway (lpoly med_growth_E_c year if year <= 1880, bwidth(10) lcolor(red)) ///
  (lpoly med_growth_E_c year if year > 1880 & year < 1920,  bwidth(10) lcolor(red)) ///
  (lpoly med_growth_E_c year if year >= 1920,  bwidth(10) lcolor(red)) ///
  (lpoly med_growth_W_c year if year <= 1880,  bwidth(10) lcolor(blue)) ///
  (lpoly med_growth_W_c year if year > 1880 & year < 1920,  bwidth(10) lcolor(blue)) ///
  (lpoly med_growth_W_c year if year > 1910,  bwidth(10) lcolor(blue)) ///
  (scatter med_growth_E_c med_growth_W_c year, mcolor(red blue)) ///
  (pci 0 1883 .6 1883, lcolor(black)) ///
  (pci 0 1918 .6 1918, lcolor(black)) if first == 1, ///
  legend(label(1 "Eastern side of TZ") label(3 "Western side of TZ") ///
  order(1 3)) ytitle("Pop. growth") 
graph export "graphs/historical_sorting/median_diff_both_curr.pdf", as(pdf) replace
restore


graph combine "graphs/historical_sorting/median_diff_1883b.gph" "graphs/historical_sorting/median_diff_1918.gph", ysize(4) xsize(8)
graph export "graphs/historical_sorting/median_growth_diff_1883_1918.pdf", replace
